arXiv:l 506.00947V 1 [stat.AP] 2Jun2015 


Localization of epileptic seizure with an approach 
based on the PSD with an autoregressive model 

ALI S. DABYEt MAHAMAT ALI ISSAKA t LAMINE GUEYE t 

November 16, 2016 


Abstract 

In this study, we present a criterion based on the analysis of EEG signals 
through the mean of the conventional power spectral density (PSD) in the 
aim to localize and detect the epileptic area of the brain. Firstly, as the 
EEG signals are commonly non stationary in practice, we processed the 
data with technique of differentiation in order to have the stationary which 
is convenient to model with autoregressive model (AR). For this, we have 
used many techniques for to determine the order which model better the 
data in this work. Therefore, we can characterize normal and abnormal 
activity which correspond to epileptic discharge for the patient. 

Our contribution in this work is the automatic detection of epilepsy seizure 
with the PSD novel approach by a better resolution in the frequency domain 
as the examination of EEG signals is often done with visual inspection of 
the rhythm (delta, theta, alpha, beta, gamma) by neurologists practitioners. 
The accuracy of the detection is estimated to 70% with the sensitivity of 
80.55% compared with the interpretation of neurologist. 
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1 Introduction 


The history of the development of time-series analysis and particularly 
the AR modeling took place in the time domain, and it began with the two 
articles of Yule |25) and Slutsky [20) in terms of theory and methodology. 
After those two articles, many others authors have worked on this area of 
mathematical statistics. These two pioneers methods were subsequently 
developed and applied by Walker (1931), Barlett (1948), Parzen (1957), 
Blackman and Tukey (1958) and Burg (1967). 

The literature on AR modeling and power spectral estimation is rather huge. 
Box and Jenkins (1970) is the first book systematically dealing with time 
series analysis within the ARMA framework. The book by Anderson (1971) 
has been written specihcally to appeal to mathematical statisticians trained 
in the more classical parts of statistics. Brillinger (1981) and Priestley (1981) 
offer wide coverage as well as in-depth accounts of the spectral analysis of 
time series. We can also cite the book of Fan and Yao (2003) who treat 
about all the area of time series with applications on simulated data. 

As journal article, according to the works who have treated the AR model 
in EEG signal, we have the hrst remarkable work of Wright J., Kydd R. 
and Sergejew A. [24| who considers the properties of parameters (natural 
frequencies and damping coefficients) obtained from segment-by-segment 
autoregression analysis of ECoG of rat. As recent works, we can cite Wang 
T.W., Guohui W. and Huanqing E. [23]. According to the authors who have 
worked on AR power spectral estimation in EEG signals, we have Abdul- 
hamit Subasi [T] who have proposed a comparative study based on the AR 
methods and EEG power spectral estimation in order to analyze and char¬ 
acterize epileptic discharges in the form of 3-Hz spike and wave complexes 
in patients with absence seizures. Also in the same main, we have O. Eaust, 
R.U. Acharya, A.R. Allen and C.M.Lin [T3] who proposed a study who deals 
with a comparative study of the PSD obtained from normal, epileptic and 
alcoholic EEG signals. The power density spectral were calculated using fast 
Eourier transform (EET) by Welch’s method, autoregressive (AR) method 
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by Yule-Walker and Burg’s method. 

In this work, the innovation relatively to [T] and m is that as the EEG sig¬ 
nals are commonly non stationary in practice, we processed the data with 
technique of differentiation in order to have the stationary which is con¬ 
venient to model with autoregressive model (AR). After this one, we used 
the AR model in order to obtain a convenient order for the model by the 
criterion of AIC, BIG and AICc. Therefore, we used the AR parametrics 
methods such as MLE, Burg’s method and Yule-Walker method in the aim 
to obtain the convenient order. 

Finally, we used our proposed methods based on the mean of the PSD to 
detect the rhythm in the frequency domain in order to discriminate epileptic 
EEG signals and normal EEG signals. 

Moreover, this approach is motivated as the analysis of EEG recording by 
neurologists practitioners is based on the visual analysis of the signals ac¬ 
cording to the rhythm delta, theta, alpha, gamma and beta which the clinical 
and physiological interests is between 0.5 and 30 Hz. This range is divided 
into four frequency bands as follows: delta (< 4 Hz), theta (4 — 8 Hz), 
alpha (8 — 14 Hz) and beta (14 — 30 Hz). 


2 Materials 

2.1 Data acquisition 

Epilepsy is one chronic and complex neurological disorder which is char¬ 
acterized by severe, possibly life-threatening episodes that can last from 
about seconds till one minute. Almost 1% of the population of the world 
suffers from the epilepsy and most are not taken care on the medical aspect. 
Besides, the unique way of highlighting the epileptic activity is the exam¬ 
ination of electroencephalography ( EEG ). The recording of the electric 
activity of the brain by the EEG is collected thanks to small electrodes 
placed on the scalp which is putted on the head of patient. 
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Epileptic seizure is an abnormality in EEG recordings and characterized by 
brief and episodic neuronal synchronous discharges with dramatically in¬ 
creased amplitude. 

Indeed, The EEG is a mean of investigation in brain diseases mainly thanks 
to its high temporal resolution allowing to detect in real time the events tak¬ 
ing place in the scale of milliseconds. It is used in cases of seizure disorders 
such as epilepsy or to report the presence of brain tumors. 



Figure 1: Original signal for normal activity 


In this work, EEG recordings of normal and epileptic subjects were ob¬ 
tained from Server of Micromed Machine in the Gentre Hospitalier Univer- 
sitaire de Fann (GHU) at Dakar. 

All the EEG signals used for this study were recorded with 128 channel 
system with 12 bit A/D resolution. In other sense, we have for each 20 sec¬ 
onds 2560 observations generated by the Micromed machine. The empirical 
analysis is based on each sequence of n = 2560 observations and the time of 
recording is between twenty minutes to one hour. 

The EEG data were obtained from two different sources, one source has 
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provided the data for epileptic seizure analysis and the other one for nor¬ 
mal subjects EEG analysis. However, we have ten patients with epileptic 
discharge and fonr which present a normal EEG signal. 



Figure 2; Original signal of patient suffering for centro-parietal and frontal 
epileptic discharge 
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3 Methods 


Let x{l),x{2), ...,x{N) the observed data of the patient which are gener¬ 
ated as times series with frequency of 128 Hz in the original EEG. At first, 
we have processed the data as stationary time series as we use the autore¬ 
gressive model (AR) which is often suitable for the condition of stationary 
in time series. For this, we use the techniques used in the literature to trans¬ 
form a non-stationary time series to stationary time series like the moving 
average, the differentiation. 

y(n) = x{n) — x(n — 1) (3-1) 

This preliminary data processing can be motivated by the fact that the EEG 
signals of the patient who present seizure are often non stationary. 


3.1 AR spectral estimation 

Given a sequence of data (a:(l), ...,x{N)) V n = 1, ...N, an autoregressive 
process of order p (abbreviated AR(p)) is a linear combination of the p most 
recent past values of itself plus an ’innovation’ term e{n) that incorporates 
everything new in the series at time n that is not explained by the past 
values. Thus, for every n, we assume that e{n) is independent of x{n — 
l),x(n-2),... 

Otherwise, a pth— order autoregressive process x(n) satisfies the equation: 

p 

x{n) =''^^a{i)x{n — i) + e{n) (3-2) 

i=l 

with e(n) ~ N{0,a‘^),a{i) G M, a(0) = 1 and a{i) G M+ are the AR coeffi¬ 
cients. 

However, our AR model can be characterized by the parameters 9 = 
{a(l), a(2),..., a(p), which we attempt to estimate by AR methods of 
estimation that we cite them in the next section. 
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Otherwise, as our goal in this study is the estimation of power spectral 
estimation (PSD), we define firstly the conventional formula of the PSD 
which is given by : 


P{f) = 


crt 


(3.3) 


\MfW 

where A{f) = 1 + a(l)e“'^^’^-^ + ... + a{p)e~P'^^^ is the hlter frequency 
response and a‘1 is the variance of the innovation term. 

We proposed an approach based on the mean of the PSD as defined by 


(3.3) which can detect the difference in the frequency resolution in order 


to classify the EEG signal according the rhythm as used by neurologists 
practitioners for interpretation. However, the proposed criterion is defined 


as 


I P{f)iiP{f)>k*P{F) 
Pm[J) = < 

I 0 otherwise 

where 


3.1.1 Maximum likelihood estimation method 

We considered that the distribution form of our given data generated 
by the Micromed machine are iid and the probability density function is 
Gaussian after confirmation with test of normality. So x ~ N{0,C{9)) and 
the PDF is expressed as 


dix]6) = - - -—-exp[— x^C ^(6*)x] 

^ ^ (27r)^/2deti/2(^(0)) 2 ^ ^ ^ 

By taking the log of (13.4^, we obtain the log-likelihood function 


(3.4) 


N N 

In d{x\ &) = -In 27r-/ 

2 2 y_i/2 


In 


cr; 






Af) 

|A(/)P 


df 


(3.5) 
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where /(/) = ;^| a:(n) exp(—j27r/n)| is the periodogram of ob¬ 

served data. 

By considering the fact that A(f) is minimum-phase, we obtain 

/■ 1/2 

/ ln\A{f)\^df = 0 
J-l/2 

and therefore 


N N N n 

Indix; 0) = -- ln27r - - - ^ y \A{f)\^df 

After differentiating and calculations, we obtain 

ri/2 


(3.6) 




1 - 1/2 


mriim 


(3.7) 


Also, after substituting and minimization, we obtain the estimated autocor¬ 
relation function defined as 


f{l) = 


wI2n=o xin)x{n+\l\), if |/|<A^-1 


0, if |/| > N 

So the set of equations to be solved in order to estimate the AR parameters 


are 


p 

d{i)f{l -i) = -r{l), I = 1, ...,p (3.8) 

i=l 

or in matrix form 


r(0) f(l) 

r(l) f(0) 

. r{p — 1) 

. f{p — 2) 


a(l) 

a(2) 


r(l) 

r(2) 

r{p — 1) f{p — 2) . 

f(0) 


_d{p)_ 


r{p) 
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(3.9) 






Thus, we obtain the so-called Yule-Walker equations and therefore we 
can estimate the AR parameters by solving this equation using Levinson 
recursion m- 


Therefore, by using the equations (3.7) and (3.8) we obtain an explicit 
form for the expressed as 

p 


al = r(0) 


(3.10) 


i=l 


Finally, we express the PSD by using those estimate parameters 9 


Pif) = 


(t: 


3.1.2 Yule-Walker’s method 


(3.11) 


The Yule-Walker AR method of spectral estimation fit the AR param¬ 
eters represented in 9 by forming a biased estimate of the signal’s autocor¬ 
relation function and a minimization of a prediction error. So, it’s directly 


based on the equation (3.9). 


However, the prediction error is expressed by 

+ 00 p 

> Lrln,'! -I- 

N 


+00 P 

T = ^ k(n)-h ^ap(i)( 


n — i) 


(3.12) 


i=l 


Therefore, by considering the equations given in (3.8) and (3.9) , the AR 


coefficients can be obtained by solving the above set of p -|-1 linear equations 
(for instance, by using the fast Levinson-Durbin algorithm) and we get 


R„a + r„ = 0 


(3.13) 


where a is the vector of AR coefficients. 

Thus, we get 

a = -Rp^Tp (3.14) 
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Once, we have the estimated AR coefficients, we can obtained the vari¬ 
ance (Tg also by using (3.8) that 


= f(0) + ^a(i)r(-i) 


(3.15) 


i=l 


Indeed, is selected in order to minimize the prediction error r 
Finally, we can expressed the PSD estimation as pQ 


PU) = 


crt 




(3.16) 


3.1.3 Burg’s method 


The Burg estimator is an innovation of the Yule-Walker estimator in the 
estimation of AR parameters. It is based on a recursive algorithm which 
is aimed at finding the sequence of values which constitute the (empirical) 
partial autocorrelation function and which are also described as reflection 
coefficients. Successive stages of the algorithm correspond to autoregressive 
models of increasing orders. At each stage, the autoregressive parameters 
may be obtained from the reflection coefficients and from the autoregressive 
parameters generated in the previous stage. For the solving of the system 
of equations, the Burg’s method used also the Durbin-Levinson algorithm 
which is the means of generating the Yule-Walker estimates recursively. 
However, the forward and backward prediction errors for a pth-order 

model are defined respectively as 

p 

£f,p{n) = x{n)+ '^ap^ix{n - i), n = p+l,...,N (3.17) 

i=l 

P 

£b,p{n) = x{n - p)+ '^a*pix{n - p + i), n = p+l,...,N (3.18) 

i=l 

where and aY = cip^p-i Vi = 1, ...,p are the respectively the forward and 
backward prediction coefficients. 


As the Burg’s method is based on minimizing the sum of the squared 
forward and backward prediction errors, we can express performance index 
(US]) by 
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(3.19) 


N 

^p= +4p(”)] 

n=p+l 

So, as the aim in Burg’s method is the estimation of reflection called for¬ 
ward and backward prediction errors, we used the Levinson-Durbin method 
and we get 


Qip^i — -|- kpQ,p—i^p—i, n — p -|- 1,..., N 


(3.20) 


^p,p—i — hp—i,p—i kpOjp—i p, n — p 1 ,..., N 

where kp is the reflection coefficient for order p. 


(3.21) 


Combining the relationship in equations (3.19), (3.20) and (3.21) leads to 


the lattice structure for computation of the forward and backward prediction 
errors, where the two prediction error series are interrelated recursively as 


£f,pin) = ef^p-i{n) + kpeb,p-i{n - 1), n = p + l,...,N (3.22) 


Sb,p{n) = eb,p-iin) + kpSf^p-i{n), n = p + l,...,N (3.23) 


Finally, the reflection coefficient kp may be chosen so as to minimize the 


equation given in (3.19), that is by setting 


dipp 

dk„ 


N 


2 E 

n=p+l 


. ^^^9eb,pin) 

+ ^b,p{n> 


dk-n 


= 0 


(3.24) 


So, partial differentiation of equations (3.22) and (3.23) with respect to h 


and substituting the results in equation (3.24), we get 


N 


^ [e/,p(n)efc,p_i(n - 1) -t- e6,p(n)e/,p_i(n)] = 0 
n=p+l 


(3.25) 
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Now substituting equations (3.22) and (3.23) in equation (3.25), we get 


N 

^ [{e/,p_i(n) + kpeb,p-iin - 1)} eb,p-i{n - 1) + {eb,p-i{n - 1) + /cpe/,p_i(n)} e/,p_i(n)] =(®.26) 

n=p-\-l 

Finally, the reflection coefficients kp can then be calculated as 


X]n=p+1 


kp = -2 


2-^n= 


n=p+l 




(3.27) 


Finally, from the estimation of AR parameters and the reflection coefficients, 
the PSD can be expressed as 


Hf) 


£ 


P 




where Sp = £f,p{n) + £b,p{n) 


(3.28) 


3.2 Optimal order determination 

The selection of order is one of the important aspects of AR method for 
power spectral density estimation. That’s because it determine the number 
of the past values requisite to predict the actual value of the time series. 
Various researchers have worked on this problem BUM and many appli¬ 
cations have been done mm- 

In this work, we used the Akaike Information Criterion (AIC) and the 
Bayesian Information Criterion (BIC) for the prediction of convenient order 
in the detection of the seizure by the PSD. In the expressions of the crite¬ 
rion used below, all the estimators are derived from the maximum likelihood 
method or its asymptotic equivalents. 


3.2.1 Akaike Information Criterion 

Akaike’s information criterion (AIC) (Akaike 1973, 1974) has been re¬ 
garded as one of the important breakthroughs in statistics in the twentieth 
century [TO] . 
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It is used to select the optimum parametric model based on observed data. 
In the context of fitting an AR model to time series data, if we regard the 


Gaussian likelihood of (3.4) as the true likelihood function, the AIC is of 


the form (after discarding some constants) 

AIC(p) = log(<T2) + ” + 

This leads to a modified form 


n 


(3.29) 


AICc(p) = log((T^) H- (3 30 ) 

n — p — 2 

The main idea is that we select the order p that minimizes AIC(p) or 
AICC(p) defined above. 

In view of the fact that the AIC tends to overestimate the orders (Akaike 
1970, Jones 1975; Shibata 1980), AICc places a heavier penalty for large 
values of p to counteract the over-fitting tendency of the AIC [10] . 


3.2.2 Bayesian Information Criterion 

Since the AIC (also AICc) does not lead to a consistent order selection 
(Akaike 1970; Shibata 1980; Woodroofe 1982), various procedures have been 
proposed to modify the criterion in order to obtain consistent estimators. 


In the same view as (3.29) and (3.31), we may define the BIC for fitting 
AR models as 


BIC(p) = log(<T2) + (3.31) 

n 
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4 Results 


4.1 Selection of optimal order 

According to the order of the model using the AIC and the BIC criterion 
with AR methods cited above represented on the Figure , we can select 
the order p = 10. 

We can remark that for the all EEG signals of epileptic patient used, the 
order became stable after p = 10. So, it’s can seen as the convenient order 
of the AR methods. 




Figure 3: Optimal order determination 
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4.2 Real EEG Data 


We have applied our algorithm of detection of severe, possibly life- 
threatening episodes of seizure in all the original EEG data. 



Figure 4: EEG signal of an unhealthy subject (epileptic patient). 



Figure 5: EEG signal of an unhealthy subject (epileptic patient). 
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I of Sample, 


Figure 6: EEG signal of a healthy subject 


4.2.1 PSD of unhealthy patients with AR methods 


In this party, we considered the EEG signal of the patient suffering 
from centro-temporal and right anterior temporal epilepsy as interpreted by 
neurologist doctors in the GHU of Dakar. 



Eigure 7: Detection of seizure by PSD according to the rhythm (red line = 
theta rhythm ,blue line= alpha rhythm). 
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Figure 8: Detection of seizure by PSD according to the rhythm (red line = 
theta rhythm ,blue line= alpha rhythm). 

In the Figure and power spectral density of an EEG signals taken 
from an unhealthy patient suffering from centro-temporal epilepsy are given. 
The detection of seizure is done according to the rhythm with the proposed 
approach. Visual inspection of this PSD conclude that we seen a delta and 
an alpha rhythm in all the derivations corresponding to centro-temporal 
seizure of epilepsy. This seem perfect according to the clinical and physio¬ 
logical interests of the rhythm as mentioned in the introduction of this study. 
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Figure 9: Detection of seizure by PSD according to the rhythm (red line = 
theta rhythm ,blue line= alpha rhythm). 
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Figure 10: Detection of seizure by PSD according to the rhythm (red line 
= theta rhythm ,blue line= alpha rhythm). 

Also, the Figure!^ and |10| show the PSD of an EEG signals of unhealthy 
patient suffering from a right anterior temporal epilepsy. We remark by 
visual examination that we seen low frequency range between theta and 
alpha rhythm on all the right temporal derivations. 
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4.2.2 PSD of healthy patients with AR methods 



Figure 11: Detection of seizure by PSD according to the rhythm (red line 
= theta rhythm ,blue line= alpha rhythm). 




Fraousncy 


Figure 12: Detection of seizure by PSD according to the rhythm (red line 
= theta rhythm ,blue line= alpha rhythm). 


Finally, the Figure and 12 show the PSD of an EEG signals of a 
healthy patient with a normal EEG. We remark by visual examination of 
those hgures that we seen null frequency range on the almost of the deriva- 
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tions and large frequency range superior to alpha rhythm. This is a coherent 
result because it’s an sleep EEG of normal patient who is an adult. 

4.3 Evaluation of the detection criterion 

In this section, the specificity, sensitivity and accuracy are used for de¬ 
termining the performance of the proposed detection method. They are 
defined as 


number of true positive decisions 

sensitivity =--- - —-— -—-, (4.32) 

number of actually positive cases 

number of true negative decisions 

specificity =--- - —-— ---, (4.33) 

number of actually negative cases 

number of correct decisions 

accuracy = ———--- - - (4.34) 

total number of cases 

We have selected three patients suffering respectively to fronto-central, 
centro-temporal and right anterior temporal seizure of epilepsy for which we 
evaluate the accuracy of the method used in this study as shown in Table 



Sensitivity (%) 

Specificity (%) 

Accuracy (%) 

Patient 1 

66.67 

63.63 

65.15 

Patient 2 

75 

50 

62.5 

Patient 3 

100 

64.24 

82.12 


Table 1: Statistical parameters classifier 
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Patient 1 

Patient 2 

Patient 3 

Fp2-F8 

0 

0 

0 

0 

1 

1 

F8-T4 

1 

1 

1 

0 

1 

1 

T4-T6 

0 

0 

1 

1 

1 

1 

T6-02 

0 

0 

1 

0 

0 

1 

Fp2-F4 

1 

0 

0 

0 

0 

0 

F4-C4 

1 

1 

1 

1 

0 

0 

C4-P4 

1 

1 

1 

1 

0 

1 

P4-02 

0 

1 

0 

1 

0 

0 

Fz-Cz 

1 

1 

1 

1 

0 

0 

Cz-Pz 

1 

1 

1 

1 

0 

1 

Fpl-F7 

0 

0 

1 

0 

0 

0 

F7-T3 

1 

0 

0 

0 

0 

0 

T3-T5 

0 

0 

1 

1 

0 

0 

T5-01 

0 

0 

1 

1 

0 

1 

Fpl-F3 

0 

0 

0 

1 

0 

0 

F3-C3 

1 

0 

0 

1 

0 

0 

C3-P3 

1 

0 

0 

1 

0 

1 

P3-01 

0 

0 

0 

1 

0 

1 


Table 2; Results of episodes of seizure detected (for each patient : first 
column= neurologist detection, second column=PSD detection) 
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AR parameters 

MLE 

Yule-Walker 

Burg 

a(l) 

0.858 

0.888 

0.858 

a(2) 

0.770 

0.713 

0.771 

a(3) 

0.049 

0.014 

0.049 

a(4) 

-0.561 

-0.524 

-0.561 

a(5) 

-0.526 

-0.451 

-0.526 

a(6) 

-0.135 

-0.117 

-0.134 

a(7) 

0.301 

0.255 

0.303 

a(8) 

0.325 

0.272 

0.325 

a(9) 

0.097 

0.099 

0.098 

a(10) 

-0.250 

-0.216 

-0.249 


0.032 

0.035 

0.031 


Table 3: AR parameters estimation of an unhealthy patient (fronto-central 
epilepsy) 


AR parameters 

MLE 

Yule-Walker 

Burg 

a(l) 

0.275 

0.269 

1.266 

a(2) 

-0.4464 

-0.433 

-0.714 

a(3) 

0.119 

0.099 

0.547 

a(4) 

-0.127 

-0.102 

-0.224 

a(5) 

0.123 

0.098 

0.224 

a(6) 

-0.144 

-0.120 

-0.243 

a(7) 

0.134 

0.116 

0.263 

a(8) 

-0.032 

-0.021 

-0.164 

a(9) 

-0.062 

-0.066 

-0.026 

a(10) 

0.108 

0.107 

0.164 


4.324 

4.372 

4.341 


Table 4; AR parameters estimation of healthy patient 
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5 Conclusion and Discussion 


In this paper, we have proposed a novel approach based on the power 
spectral density of the EEG signals according to the range of the frequency. 
We used the AR methods for PSD such as Maximum likelihood Estimator, 
Burg’s method and Yule-walker’s method in order to estimate the param¬ 
eters of the PSD. Moreover, we used the Akaike Information Criterion 
(AIC) and the Bayesian Information Criterion (BIC) for the prediction of 
convenient order in the detection of the seizure by the PSD. 

However, we have presented the PSD of EEC signals taken from an 
unhealthy patient and a healthy for the detection of the seizure. By Visual 
inspection of those PSD, we remark that we have a low frequency around 
the rhythm alpha for epileptic patient, and a null frequency and large 
frequency for the healthy patient used in our study. This confirm the 
clinical and physiological range frequency used by neurologists doctors for 
visual inspection of EEC signals. 

This contribution is very helpful for neurologists practitioners because it 
improve the detection in frequency domain giving a better resolution in the 
sense of reduction of time detection as the examination of EEC signals is 
often done with visual inspection of the rhythm. 

As prospects for the future, it will be interesting to get a model who permit 
as to detect the seizure according to the change in frequency and temporal 
domain with non parametric methods. 
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